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We perform the first fully nonlinear numerical simulations of black-hole binaries with mass ratios 
100:1. Our technique is based on the moving puncture formalism with a new gauge condition and 
an optimal choice of the mesh refinement. The evolutions start with a small nonspinning black hole 
just outside the ISCO that orbits twice before plunging. We compute the gravitational radiation, 
as well as the final remnant parameters, and find close agreement with perturbative estimates. We 
briefly discuss the relevance of these simulations for Advanced LIGO, third-generation ground-based 
detectors, LISA observations, and self-force computations. 
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Introduction: The orbital evolution and computation of 
gravitational radiation from black-hole binaries (BHB) in 
the small-mass-ratio limit remains one of the most chal- 
lenging problems in General Relativity. This was rec- 
ognized early on by Regge and Wheeler over 50 years 
ago pQ. Zerilli then completed the formulation of the 
first order perturbations around a Schwarzschild BH in 
1970 [2j. Three years later, Teukolsky [3] provided a new 
formalism to study perturbations around Kerr BHs. In 
order to take into account the decay of the orbit of the 
small BH due to the emission of gravitational radiation, 
second order effects have to be included in those compu- 
tations. This problem turned out to be very challenging, 
and only since 1996 [31 [5] has there been a consistent for- 
malism for the "self-force" corrections to the background 
geodesic motion of a small BH orbiting around a larger 
one. The explicit implementation of such formalism into 
a computational scheme remains challenging, although 
recent progress along this line is encouraging [6]. 

The dramatic breakthroughs in the numerical tech- 
niques to evolve BHBs 7HH] transformed the field of Nu- 
merical Relativity (NR) and we are now in a position 
to evolve binary systems in the intermediate mass ratio 
regime. Two years ago the merger of spinning |10] bina- 
ries with mass ratio q = toi/to2 = 1/8 and nonspinning 
binaries [11 j with q = 1/10 were published. More recently, 
detailed long term evolutions of BHBs with q = 1/10 
and q = 1/15 were studied and validated against per- 
turbation theory (T2l [13] . In this Letter we present the 
first fully nonlinear numerical simulations of the merger 
of small-mass-ratio BHBs. As a case study, we evolve a 
nonspinning BHB with mass ratio q = 1/100 for over two 
orbits prior to merger, and resolve the entire waveform 
for three grid resolutions, proving numerical convergence 
of the results. The success of our approach is based on 
enhancements of the moving puncture numerical tech- 
niques that adapt the gauge and grid structure to the 



TABLE I: Initial data parameters. The punctures are located 
on the x-axis at positions x\ and X2, with puncture mass 
parameters (not horizons masses) mi and rri2, and momentum 
±p. The punctures have zero spin. The ADM mass Madm is 
1M and q = 0.01000004. 



Xl 


4.95256 


X2 


-0.0474374 


Px 


-0.0000102652 


Py 


0.00672262 


mi 


0.00868947 


Vtl2 


0.989619 



small-mass-ratio limit. 

The techniques described in this Letter can be used 
in the spinning BHB case and for even smaller mass ra- 
tio inspirals. This has important consequences for as- 
trophysics and gravitational wave observatories such as 
the second-generation Advanced LIGO detector, third- 
generation ground-based detectors, and LISA. Supermas- 
sive BH collision at cosmological scales are most likely to 
occur in the mass ratio range 1:10 - 1:100 14] and will be 
observable by LISA, while collision of intermediate mass 
BHs and solar mass BHs will lie in the sensitivity band of 
second and third generation ground-based detectors [TBT - 

E]. 

Fully Nonlinear Numerical Simulations: In Table [I] we 
give the initial data parameters for our q = 1/100 
BHB simulations. We evolved this BHB data set us- 
ing the LazEv |18, implementation of the moving punc- 
ture approach [8, 9J. Our code used the Cactus/Einstein 
toolkit [HI [20] and the Carpet [21] mesh refinement 
driver to provide a 'moving boxes' style mesh refinement. 
We use AHFinderDirect [3J] to locate apparent hori- 
zons. We measure the magnitude of the horizon spin 
using the Isolated Horizon algorithm detailed in [55] . 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with a 
modified 1+log lapse and a modified Gamma-driver shift 
condition 0H3], and an initial lapse a(t = 0) = 2/(1 + 
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ip% L ). The lapse and shift are evolved with (dt—^d^a = 
-2aK, d t (i a = |f a - r](x k ,t)/3 a , where different func- 
tional dependences for rj(x k ,t) have been proposed in 
[T51 I25T - I2!?] . Here we use a modification of the form pro- 
posed in [25], r](x k ,t) = Roy/fidiWdjW/Hl - W a ) b ), 
where we chose i?o = 1.31 and W is the evolved con- 
formal factor. The above gauge condition is inspired 
by, but differs from Ref. [25J between the BHs and in 
the outer zones when a ^ 1 and 6 / 2. Once the 
conformal factor settles down to its asymptotic i/> = 
C '/ \fr + O(l) form near the puncture, r\ will have the 
form r\ — (R /C 2 )(l + b(r/C 2 ) a ) near the puncture and 
r\ = Ror b ~ 2 M I (aM) b as r — > oo. Our exploration of 
the (a, b) parameters showed that the (1,2) case leads 
to numerical instabilities on coarse grids, while the (2,1) 
and (1,1) cases lead to noisy waveforms and slower gauge 
speeds. In practice we used a = 2 and b = 2, which 
reduces 77 by a factor of 4 at infinity when compared to 
the original version of this gauge proposed by [25]. We 
note that if we set b = 1 then 77 will have a 1/r falloff at 
r = 00 as suggested by [28] . 

In order to chose the width of the refinement levels clos- 
est to the small BH, we examine the potentials for pertur- 
bations about a nonspinning BH. The idea is that we need 
to model the curvature and the gravitational radiation 
emitted by the small BH (which drives the merger, and 
hence the physics) . At the zeros of the derivative of the 
potentials, the variations are minimized. Furthermore 
the separations between zeros increases, naturally lead- 
ing to a choice of small-width, high resolution grids be- 
tween the first zeros, one step lower in resolution between 
the second two, followed by a sequence of coarser grids. 
According to Chandrasekhar [30] pl60] the even/odd (±) 
parity effective potentials of a Schwarzschild BH can be 
written as 



V? 



±6M 



4f_ 

dr* 



+ (6M) *f* + 4A(A + 1) /, (1) 



where 



/ = 



(r - 2M) 
2r 2 (Ar + 3M)' 



A = 



+ 2)(£-l). (2) 



Note that both potentials are numerically very close to 
each other, hence we consider the vanishing of the deriva- 
tive of the average of the two (in isotropic coordinates R) 
when constructing the grid. For 1 = 2 and M = 1 this 
takes the explicit form 



d(V^ 



v- 



-384 
(64i? 6 



®R f=2, M=l 

R (2 R - 1) (16 i? 4 - 4 R 3 - 60 R 2 - R + l) 



(3) 



(2R+ l) 9 (4i? 2 + 10i? + l) 3 
288i? 5 + 480i? 4 + 256i? 3 + 120i? 2 + 18 R + l) 



Ideally, we would like to place the AMR boundaries 
around the small BH near the zeros of this function 



(R/m = 0.0, 0.1207998431, 0.5, 2.069539112) The loca- 
tion of the zeros for £ > 2 changes little from the above 
figures. Since we do not want to over-resolve the inte- 
rior, this suggests that the first grid level should cover 
the whole small hole up to its horizon, and the next grid 
level up to 4 times the initial horizon radius of the small 
hole in the initial quasi-isotropic coordinates. 

In Ref. |13j we provide an alternative method of ex- 
trapolation of waveforms based on a perturbative prop- 
agation of the asymptotic form of ip4 at large distances, 
leading to the following simple expression 



lim [rVl m (r,t)] 



(4) 







+ 0(rol), 



where robs is the approximate areal radius of the sphere. 
This formula is applicable for robs ^ 100M. And note 
that it is also important to remove the low frequency 
components [31] in ^4 (since it is inside an integral). 

Results and Analysis: Our simulation used 15 levels of re- 
finement (around the smaller components), with central 
resolutions as high as M/7078, and 9 levels of refinement 
around the larger component. The outer boundaries were 
located at 400M and the resolution in the boundary zone 
was h = 2.3148M for our finest resolution run. The BHB 
performs ~ 2 orbits prior to merger [as seen by the for- 
mation of a common apparent horizon (CAH)], which 
occurs roughly 160M after the start of the simulation. 
In terms of computational expense, a medium resolution 
run requires 500,000 SU and approximately one month 
of runtime. In order to reduce the total runtime, we used 
an aggressive choice of CFL (dt = l/2h), which leads to 
significant BH mass loss when compared to dt = 1/Ah. 

Table |TT] shows the results of evolution. We note that 
the smaller BH mass is conserved to within 0.23% dur- 
ing the inspiral and plunge phases, while the mass of 
the larger BH is conserved to within 0.003%. In Fig. [I] 
we show the xy projection of the orbital trajectories for 
the two highest resolution runs. From the figure we can 
see that the initial jump in the orbit pushes the binary 
slightly outside the ISCO, leading to an additional orbit. 
In Fig. [2] we show the orbital radius as a function of time 
and resolution. Note that the orbital radius supercon- 
verges at low resolution and converges quadratically at 
high resolution ( this quadratic error may be due to time 
prolongation effects, as well as effects due to an aggres- 
sive choice of CFL factor). In Fig. [3] we show amplitude, 
as well as phase convergence, of the (£ = 2, m = 2) mode 
of "04. The reduced order of convergence is due to an 
aggressive choice of CFL factor. 

The apparent superconvergence in the trajectories and 
waveforms when considering the three coarsest resolu- 
tions is indicative that the lowest resolution is just en- 
tering the convergence regime. That is, this resolution 
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TABLE II: Remnant horizon parameters and radiated energy- 
momentum. Here we provide SM^ = Madm — Mh and 
SS^j = J adm — Sh, which are small numbers obtained by 
taking the difference between two much larger numbers. The 
calculation of 8S'% is relatively inaccurate because it requires 
an extrapolation to infinite resolution. 



10 E rad 


6.0 ±0.1 


10 5 5M* H 


7.0 ± 1.0 


100a 


3.33 ±0.02 


W 4 J ra d 


5.0 ±0.2 


IQ A 5S* H 


3.0 ±2.0 


Vkick 


1.07 ± 0.05km s~ 
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FIG. 1: An xy projection of the trajectories for the two high- 
est resolutions of the q = 1/100 configuration. The dotted 
circle corresponds to the ISCO radius while the small filled- 
in circle corresponds to the point on the trajectory where a 
common horizon is first detected. Note the initial "jump" in 
radius (see Fig. [5| due to the initial data radiation content. 




|<Kh,/1 -2)-iMh/1 2 )]"1 396 (2nd) 
[♦(tyi.2)-t(h /1.2')p.16(1st) 




FIG. 3: Convergence of the amplitude and phase of the (£ = 
2,m = 2) mode of ip4. The phase converges to second order 
prior to the peak in the amplitude. The vertical line shows the 
point when uj — 0.2. Note the good agreement in amplitude 
(the curves have been translated). The phase error at cj = 0.2 
is 0.44 radians. 




cannot be far from the convergence regime because all 
four resolutions lie in a monotonic convergence sequence. 
And importantly, the deviations between the next three 
resolutions are very small compared to the deviation be- 
tween the lowest two resolutions, indicating that these 
three resolutions are safely inside the convergence regime. 




FIG. 2: The orbital radius as a function of time and resolution 
for the q — 1/100 configuration. Note the initial "jump" in 
the orbit due to the initial data. 



FIG. 4: The remnant spin a/Mn and the radiated mass from 
infinite separation 8M/M, as a function of the symmetric 
mass ratio r\ = q/(l + q) 2 , for q — 1/10, 1/15, 1/100, as well 
as the empirical formula prediction. 



Finally, in Fig. [4] and Tabic [TTT we show the remnant 
spins and total radiated mass as a function of mass ratio 
[13] for q = 1/10, 1/15, 1/100 and the predictions based 
on our empirical formula [32]. Note that no fitting is 
involved in this figure. 

The amount of energy and angular momentum ra- 



TABLE III: Remnant spin and total radiated mass (starting 
from infinite separation) as a function of mass ratio q as mea- 
sured in our simulations and as predicted by our empirical 
formulae 1321. 



q 


1/10 


1/15 


1/100 


a (Computed) 


0.2603 


0.18875 


0.0333 


a (Predicted) 


0.2618 


0.1903 


0.03358 


SM (Computed) 


0.00826 


0.00507 


0.000618 


SM (Predicted) 


0.00806 


0.00498 


0.000604 
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diated when the (2, 2) mode frequency is larger than 
M«2,2 > 0.167 is given by (adding up to i — 4 modes) 
5E/M = 0.000047 ± 0.000001 and SJ/M 2 = 0.00034 ± 
0.00001, which agrees to within 4% with the particle limit 
predictions of SE/M = 0A7n 2 and SJ/M 2 = 3.Un 2 [33]. 

Conclusions and Discussion: We have successfully 
evolved a 1:100 BHB system for the last two orbits before 
merger and down to the final Kerr hole remnant. We have 
achieved this within the moving punctures approach by 
adapting the gamma-driver shift condition with a vari- 
able damping term. Also crucial for evolutions is an op- 
timal choice of the mesh refinement structure around the 
small BH. We used the Regge- Wheeler- Zerilli potentials 
to guide the setting up of the initial grids. This helps 
optimizing the large resources required to evolve small 
q binaries. The numerical convergence of the waveforms 
displayed here, and the successful comparisons with per- 
turbative results [12l [13] , show this approach is validated 
in the intermediate mass ratio regime, and can be applied 
to even smaller g's (and larger initial separations into the 
post-Newtonian regime) . 

The feasibility of simulating extreme mass ratios by 
purely fully nonlinear numerical methods, as demon- 
strated in this work, allows us to look more optimistically 
at the task of generating a bank of templates for second 
and third generation ground-based detectors and LISA. 
Methods like those described in [TJJ[T3], that combine NR 
and perturbative techniques can be used to speed up the 
generation of those templates. And finally we now also 
have a direct way of validating self-force computations 

The techniques presented here would appear to apply 
in a straightforward manner to even smaller mass ratios 
q and to initially spinning BHs. Fine tuning of the qua- 
sicircular orbital parameters plays an important role in 
preparing these runs, given the very low level of gravi- 
tational radiation they generate. So far we see that the 
method [33] developed for equal mass BHBs to lower the 
eccentricity seems to work, but it requires extra runs for 
initial experimentation. Hence it would be important to 
evolve initial data with lower spurious radiation content 
and some true inspiral wave information 35 . 
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